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Quantized vortices stunningly illustrate the coherent nature of a superfiuid 
Bose condensate of alkali atoms. Introducing an optical lattice depletes this co- 
herence. Consequently, novel vortex physics may emerge in an experiment on a 
harmonically trapped gas in the presence of a rotating optical lattice. The most 
dramatic effects would occur in proximity to the Mott state, an interaction dom- 
inated insulator with a fixed integer number of particles per site. We model such 
a rotating gas, showing that the lattice-induced spatial profile of the superfiuid 
density drives a gross rearrangement of vortices. For example, instead of the 
uniform vortex lattices commonly seen in experiments, we find parameters for 
which the vortices all sit at a fixed distance from the center of the trap, forming 
a ring. Similarly, they can coalesce at the center, forming a giant vortex. We 
find that the properties of this system are hysteretic, even far from the Mott 
state. We explain this hysteresis in terms of vortex pinning, commensurability 
between vortex density and pinning site density, and energy barriers against 
changing the number of vortices. Finally, we model time-of-fiight expansion, 
demonstrating the experimental observability of our predictions. 

Atomic clouds in a rotating optical lattice provide a system which is at the intersection of 
several major paradigms of condensed matter physics. These rotating clouds may display a 



a 



superfiuid phase, an interaction driven insulating phase [l[, and even analogs of the fractional 
quantum Hall effect 21]. Here we explore the theory of vortices in such systems. Most 
importantly, we investigate how proximity to the Mott insulator phase impacts the vortex 
configurations. 

Considering a uniform gas of atoms of mass m in an optical lattice rotating with frequency 

Q, there will be three macroscopic length scales in the problem: the lattice spacing d, the 

magnetic length i = ^Jhjm^^ and the particle spacing n~^/^, where fi = h/2'n is Planck's 

constant. Even without interactions, the commensurability of these various lengths can 

ead to nontrivial physics - the single particle spectrum, the Hofstadter butterfly, is fractal 

3l. For a system of interacting bosons, this fractal spectrum leads to a modulation of the 

nn 

boundary between superfiuid and Mott insulating phases [^, 15|]. Further, the vortices in a 
superfiuid on a rotating lattice develop extra structure: for example they can have cores filled 
with the Mott state [6i]. Related to this structure, the vortices can become pinned at optical 
lattice potential minima, in contrast to the more usual pinning at potential maxima [4f|. We 



explore how this physics plays out in a finite system, as can be studied in experiments. 



We consider a harmonical 
lattice. A recent experiment 



y trapped gas of bosons in a rotating two-dimensional square 
71] realized exactly this scenario by placing a rotating mask in 
the Fourier plane of a laser beam which forms an optical dipole trap. The mask contained 
three/four holes, producing a triangular/square lattice in the image plane, where the atoms 
were trapped. The lattice spacing was large due to the nature of their optics but can be made 
comparable to the wavelength of the trapping light. If one reduced this lengthscale, one could 
explore the tight binding limit, where atomic motion is limited to the lowest band, and Mott 
insulating physics could become important. We consider the theory of this tight binding 



lich introduce 
9, 



10 



11|. We 



limit. This limit may also be reached through quantum optics techniques^ 

phases on the hopping matrix elements for atoms in a non-rotating lattice 

choose to study a two-dimensional cloud, as it provides the simplest setting for investigating 

vortex physics. This is also an experimentally relevant geometry, as the dimensionality of 

the system can be controlled by using an anisotropic harmonic potential, o r op tical lattice. 



where the hard trapping direction is along the optical lattice rotation axis 12] 



In the tight-binding limit, this system is described by the rotating Bose- Hubbard hamil- 



tonian 



13|: 



H = -^ (tijalcij + h.c.j + ^ ( —hi {hi - 1) - ^lihi \ (1) 



where tij = texp i JI' df- A{r) is the hopping matrix element from site j to site i. The 
rotation vector potential, which gives rise to the Coriolis effect, is A(r) = [m/K) iVtxr] = 
TTu [xy — yx), where z/ is the number of circulation quanta per optical-lattice site. The 
local chemical potential fii = fiQ — m {u^ — Q"^) rf/2 includes the centripetal potential. In 
these expressions, /xq is the central chemical potential, u is the trapping frequency, fi is 
the rotation speed, Tj is the position of site i, m is the atomic mass, a\ (di) is a bosonic 
creation (annihilation) operator, hi = a\di is the particle number operator for site i, and U 
is the particle-particle interaction strength. The connection between these parameters and 
the laser intensities are given by Jaksch et al. [l3|. Here, and in the rest of the paper, we 
use units where the lattice spacing is unity. 

As described in the Methods section, we use the variational Gutzwiller ansatz to calculate 
the properties of this system. This approach describes well both the superfluid and Mott 
insulating states, and has the important property that the superfiuid density p^ is not equal 



to the total density p. This is one of the key ways in which the lattice system differs from 
a traditional weakly interacting gas of bosons, and much of the physics we see is contingent 
upon this feature. 

We systematically explore the phase space, varying the parameters in the hamiltonian. 
We simulate clouds with diameter from 30-60 sites, comparable to the sizes studied in 
experiments [ij, [l5|. For the largest simulations we impose four-fold rotational symmetry, 
but introduce no constraints in the smaller simulations. We find interesting results both 
near and far from the Mott regime. For weak interactions the dominant physics involves 
the pinning of vortices between lattice sites, while near the Mott boundary, the non-uniform 
superfluid density becomes central. 

Commensurability, Hysteresis, and Pinning 

Previous work, focusing on the multi-band weak lattice limit, found results similar to 



those we see in our tight binding model far from the Mott regime la, llTJ]. We find a great 
variety of vortex structures, including patterns resembling square vortex lattices. These are 
most stable at the rotation rates where they are commensurate with the underlying optical 
lattice. Commensurate Bravais lattices exist when 1/z/ is an integer, and commensurate 

n n n 

square lattices when v = Xjir? + w?\ for integral n and m [7|, ll6|, ll7|]. Which vortex 
patterns appear in a simulation, or in an experiment [18|, depends on how the system is 
prepared. This hysteresis occurs because the energy landscape has many deep gorges with 
near-degenerate energies, separated by high barriers. 

To illustrate this energy landscape, we fix tjU = 0.2 and po/U = 0.3, and study how the 
energy evolves as we change the rotation speed. First, starting with the non-rotating ground 
state, we sequentially increase the rotation speed, using the previous wavefunction as a seed 
for our iterative algorithm. We adjust u as we increase Q so that the cloud size, related 
to the Thomas- Fermi radius, Rtf = \/— r^zo^, remains effectively fixed. The energy as a 



function of rotation speed, shown in Figure [T]^a), has a series of sharp drops, corresponding 
to the entry of one or more vortices from outside of the cloud. At these rotation speeds the 
system jumps from one local minimum of the energy landscape to another. 

Figure [T] (b)-(g) shows the superfluid density and phases associated with the vortex 
patterns found during this adiabatic increase in rotation speed, where we impose four-fold 
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FIG. 1: Adiabatic spin-up (color, one-column). Properties of cloud during adiabatic spinup with 
parameters [tjU = 0.2, fj,o/U = 0.3, Rtf = 15). (a) Energy vs rotation rate. Sharp drops indicate 
vortex formation. Energy scaled by on-site interaction parameter U, and rotation rate quoted as 
z^, the number of circulation quanta going around each plaquette. (b), (d), (f) Superfluid density 
profile at parameters labeled in (a). Dark-to-light shading corresponds to low-to- high density, and 
position is measured in lattice spacing. Dark spots correspond to vortex cores. Red and green lines 
are guides to the eye. (c), (e), (g) Superfluid phase represented by Hue. Solid white circle denotes 
edge of cloud. Dashed lines are guides to the eye. Black circles denote vortex locations. In (b), 
(c) and (f), (g) rotation speed should favor a commensurate square vortex lattice rotated by 7r/4 
from the optical lattice directions, (d), (e) represents an incommensurate rotation speed. 



rotation symmetry about the trap center. Subfigures (b) and (c) show a regular square 
vortex lattice seen near the commensurate v = 1/(2 x 6^). Subfigures (d) and (e) show the 
vortex configuration at z/ ~ 1/(2 x 3.76^) which is intermediate between the commensurate 
values V = 1/(2 x 3^) and u = 1/(2 x 4^). Rather than forming a square pattern, the vortex 
configuration appears to be made of concentric rings. Such ring-like structures also occur 
for superfiuids rotating in hard- walled cylindrical containers p^], where boundaries play an 
important role. Changing the shape of the boundary should change the vortex configuration. 
For example, by using an asymmetric trap one should produce elliptical shells of vortices. As 
one increases u towards u ~ 1/(2 x 3^), a domain containing a square vortex lattice begins 
to grow in the center of the trap. As seen in subfigures (f) and (g), at commensurability the 
domain only occupies part of the cloud, even though one would expect that a uniform square 
lattice would be energetically favorable. The inability of the system to find the expected 
lowest energy configuration during an adiabatic spin-up is indicative of the complicated 
energy landscape. 

The patterns which we find are largely determined by the symmetry of the instabilities 
by which vortices enter the system. For example, even when we do not impose a four-fold 
symmetry constraint this adiabatic spin-up approach never produces square vortex lattices 
oriented at an angle other than tt/4 with respect to the optical lattice axes. On the other 
hand, we readily produce other commensurate vortex lattices by choosing the appropriate 
rotation speed and seeding our iterative algorithm with the corresponding phase pattern. 
We have verified this approach with square vortex lattices oriented at various angles with 
respect to the optical lattice, taking z/^^/^ = {4, 5, 6, 7, 8, 9, 10}, (5i^)"^/^ = {2, 3, 4, 5, 6} and 
(10z/)-V2 = 12,3,4}. 

We further explore the history dependance of the vortex configurations by increasing, 
then decreasing Q. We do not impose a four-fold symmetry constraint, but take a smaller 
system with Rtf = 7. At any given Q, the energy shown in figure [2] (a) depends on the sys- 
tem's history. The increasing(blue)/decreasing(red) rotation curve has sharp energy drops 
signaling the introduction/ejection of vortices to/from the system. The energy drops oc- 
cur at different Q for spin-up and spin-down, indicating that the critical rotation speed 
for a vortex to enter or exit the system is different. Generically, there are more vortices 
in the system on spin-down than on spin-up. Depending on Q, one may find a lower en- 
ergy state by increasing (subfigs. (b) and (c)) or by decreasing (subfigs. (d) and (e)) the 
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FIG. 2: Hysteresis (color, one-column), (a) Energy versus rotation rate during increase (blue 
line) and decrease (red line) of v. Energy steps in the blue (red) curve correspond to nucleation 
(expulsion) of vortices, (b)-(e) Order parameter complex phase for parameters labeled in (a). Black 
circles are drawn around vortex cores, and white circles indicate the approximate extent of the gas. 

rotation rate. As demonstrated by subfigs. (d) and (e), vortex configurations produced 
during spin-up typically have the four-fold rotational symmetry of the optical lattice, while 
the vortex configurations calculated during spin-down are more likely to break this symme- 
try. An experiment will display the same qualitative features, but slightly different vortex 
configurations. 

Influence of Mott physics 



Next we investigate the system in a deeper optical lattice, close enough to the Mott 
insulator phases so that p'^ and p are noticeably different. A basic understanding of this 
region comes from investigating the phase diagram of the homogeneous system, where fi = 



uj = 0. The mean-field phase diagram for the uniform non- rotating system is shown in 
Fig. [31 Contours of fixed p and p'^, are indicated by red and black curves. The superfiuid 
density vanishes in the Mott regions where the density is commensurate with the lattice, 
and the number of particles per site labels each lobe. For a sufficiently gentle trap, the gas 
looks locally homogeneous, and its density at any point r can be approximated by that of a 
uniform system with chemical potential p{r) = po ~ V{r). We go beyond this local density 
approximation (LDA) in our calculations, however it is useful for gaining qualitative insight 
into the density profiles we expect. 

For illustration, consider the three vertical (fixed t/U) lines in figure [3l corresponding to 
the LDA trajectories of the chemical potential in three separate trapped clouds. Trajectory 
"Tl" represents a system consisting of alternating superfiuid and Mott-insulator phases. 
Trajectory "T2" represents a Mott-insulating core surrounded by a superfiuid ring. Trajec- 
tory "T3" represents a system with a high p'^ central region surrounded by a plateau of near 
constant p'^. We self-consistently calculate the density profile of atomic clouds where po/U 
corresponds to the tops of the T3 and T2 lines. 

In figure H] we display the density and complex phase profiles corresponding to T3 for 
the non-rotating and rotating optical lattice cases. These calculations are performed at 
(t/U = 0.06, po/U = 0.7, Rtf = 10), and with four-fold rotation symmetry enforced. The 
cloud has a diameter of roughly 26 sites. In the non-rotating case we see that the superfiuid 
density profile (subfigure (a)) has the plateau structure predicted by the LDA. No such 
plateau is seen in the total density (subfigure (c)). 

The phase of the superfiuid order parameter is uniform in the absence of rotation (sub- 
figure (e)). Starting from this non-rotating configuration we gradually increase the rotation 
speed to z/ = 0.04. Rather than forming a square lattice, the resulting vortices form a ring 
around the central p'^ peak (subfigure (c)). This configuration is favored because it minimizes 
the sum of competing energy costs: the rotation favors a uniform distribution of vortices, 
but the single vortex energy is smallest where p'^ is low. 

As seen in subfigure (f), the phase of the superfiuid inside the ring is essentially constant. 
This can be understood by an analogy with magnetostatics. Given that the velocity field 
around a vortex line, v, is analogous to the magnetic field around a current carrying wire, it 
obeys the equivalent of Ampere's law §w-M = {h/m)Ny, where N^ is the number of vortices 
enclosed in the contour of integration. In the limit where one can neglect the discreteness 
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FIG. 3: Mean-field 2D Bose-Hubbard model phase diagram (color, one-column). Red 
(black) curves: contours of constant total-particle (superfluid) density. Dark-to-light purple shad- 
ing: low-to-high superfluid density. Qualitative understanding of trapped gas is produced by 
following lines of fixed t/U, such as those illustrated by orange lines {Tl, T2, T3} discussed in the 
text. 

of the vortices in the ring, one finds that the fluid inside is motionless, while the fluid 
outside moves with the same azimuthal velocity that one would find if all the vortices were 
at the geometric center of the cloud. Even with only eight vortices our system appears to 
approach this limit. If one increased the rotation speed, one could imagine finding a state 
with several concentric rings of vortices in the plateau. Similarly, increasing fio/U could 
produce multiple superfluid plateaus, each of which could contain a ring of vortices. This 



A^~ - 





FIG. 4: Ring vortex configuration (color, one-column). (a), (c), (e) Superfluid 
density, total-particle density, and complex phase profile for a non-rotating system at 
(t/U = 0.06, fio/U = 0.7, Rtf = 10) corresponding to T3 in figure H (b), (d), (f) Same after 
adiabatically increasing rotation rate to ly = 0.04. The vortices sit in a circle surrounding the 
central high superfluid density region. The central fluid is approximately stationary with constant 
phase, while the outer velocity is approximately azimuthal. 



structure is reminiscent of Onsager and London's original proposal of vortex sheets in liquid 
helium |2o|. 

In figure [5] we investigate the state corresponding to trajectory T2, where the LDA 
predicts a superfluid shell surrounding a Mott core. Rather than forming a lattice of discrete 
vortices, one expects that this system will form a "giant" vortex when rotated: the vortices 
occupying the Mott region, leaving a persistent current in the superfluid shell. The energy 
barriers for changing vorticity should be particularly high, so we generate the rotating state 
in two stages. (A similar protocol could be used in an experiment.) We start with a non- 
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rotating system (z/ = 0) at weak coupling {t/U = 0.2,/io/f/ = 0.3), gradually increasing the 
rotation to z/ = 0.032, where we find the square vortex lattice illustrated in subfigures (a), 
(c) and (e). We then adiabatically reduce t/U from 0.2 to 0.03. As we reduce t/U, the 
central p^ drops, while p approaches unity there. Eventually we see a Mott regime at the 
center of the cloud. During the evolution, we find that 8 of the vortices migrate to the edge 
of the trap and then escape, while four of the vortices move towards the center. These four 
vortices coalesce at the center of the trap and effectively form a vortex of charge 4, which is 
illustrated in subfigures (b), (d) and (f). Such a dense packing of vorticity would be unstable 
in the absence of the optical lattice. For larger systems with higher rotation rates we have 
produced giant vortices with much higher quantization. 

Due to its multiply connected topology, a ring is one of the archetypical geometries used in 



theoretical discussions of superfiuidity 2l|, |22|, |23| . There are several experimental schemes 



for creating a ring-shaped trap (for example using the dipole forces from a blue detuned laser 
focussed at the center of an atomic cloud) [2J], and many theoretical studies of g iant vortex 



formation stabilized by a quadratic-plus-quartic potential 25 



26 



27 



28 



29|. Here the 



multiply connected geometry is spontaneously formed by the appearance of the Mott state 



30|, 



in the center of the cloud. As was found in a related study by Scarola and Das Sarma 
this Mott region effectively pins the vortices to the center. 

By changing t/U one may study a few other interesting structures. For example, one can 
engineer a situation where a central superfluid region is surrounded by a Mott ring followed 
by a superfluid ring. At appropriate rotation speeds one produces a conflguration which has 
properties of both the states seen in figure H] and in figure [5l One will find no vortex cores 
(all of the vorticity is confined to the Mott ring), the central region will be stationary, and 
the outer region rotates. 

Another interesting limit is found when one decreases the thickness of a Mott/superfiuid 
region so much that it break up into a number of discrete islands. Small Mott islands act 
as pinnin g ce nters, while small superfiuid islands form an analog of a Josephson junction 



array [3ll, l32 1 . 
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FIG. 5: Giant vortex (color, one-column), (a), (c), (e) Superfluid density, total-particle density 
and complex phase field for a square vortex-lattice state prepared by adiabatically increasing 
rotation rate at relatively low optical lattice depth {t/U = 0.2,fj,o/U = 0.3, Rtf = 7, u = 0.032). 
(b), (d), (f) Same after adiabatically increasing optical lattice depth so that t/U = 0.03. A 
superfluid ring surrounds a Mott state, and four of the vortices have coalesced to form a giant 
vortex of charge 4. 

Detection 

Many of the structures we have discussed would be very hard to detect using in-situ 
absorption imaging. As is exemplified by figure 4(d), the vortices do not necessarily have 
a great influence on the density of the cloud. This is principally because near the Mott 
boundary the superfluid fraction becomes small: even though the superfluid vanishes in the 
vortex core, this does not appreciably change the density. Two other pieces of physics also 
influence the visibility. First, near the Mott boundary one can produce vortices with Mott 
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cores [Gl . Depending on the bulk density, this can lead to vortices where there is no density 
suppression at all, or even a density enhancement. Second, the lengthscale of the vortex 
core, the superfluid healing length, varies with U/t. For both very large and very small 
U/t the healing length is very large, while at intermediate couplings it is comparable to the 
lattice spacing, possibly below the resolution of typical optics. 

We argue that the vortex structures will be much more easily imaged after time-of- 



flight expansion of the cloud 33|, |3J]. The density after time-of-flight expansion is made 
of two pieces. There will be a largely featureless incoherent background from the normal 
component of the gas, and a coherent contribution from the superfluid component. Due 

to interference from the different lattice sites, the coherent contribution will form a series 

ii 

of Bragg peaks [1]. Each peak will reflect the Fourier transform of the superfluid order 
parameter. In the methods section we present a simple approach, neglecting interactions 
during the time-of-flight, which allows us to calculate the time-of-flight expansion images. 
Figure [6] illustrates the density pattern which will be seen if the rotating cloud in fig. H] is 
allowed to expand. 

Summary 

In this paper we have studied the vortex configurations in a harmonically trapped Bose 
gas in the presence of a rotating deep optical lattice. Far from the Mott phases we illustrated 
the effects of vortex-pinning and commensurability in determining ground state vortex con- 
figurations. We studied the hysteretic evolution of vortex configurations as the rotation 
speed is varied. We explored the influence of the Mott phase on the rotating superfluid, 
encountering a series of novel vortex structures such as rings and giant vortices. Finally, we 
explained how these structures can be observed in time-of-flight images of the atomic cloud. 
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FIG. 6: Time-of-flight expansion (color, one-column), (a) Column density p scaled by the 
central column density pQ as a function of space for the state in fig. U] (b) , (d), (f) after expanding 
for time t. Positions are measured in terms of scaling parameter Dt = ht/mX, where A is the size 
of the initial Wannier state in each site of the lattice. This scaling solution holds in the long time 
limit where Dt is large compared to the initial size of the cloud, (b) One dimensional cut through 
center of (a). Note the incoherent background between the Bragg peaks, (c) Close-up of the central 
Bragg peak, corresponding to the Fourier transform of the superfluid order parameter. There are 
8 dips in the outer crest, corresponding to the 8 vortices in the initial state. 

Methods 

General 

Both the superfluid and Mott insulator are well described by a spatially inhomogeneous 
Gutzwiller product ansatz [13i] . 



M 



\'^Gw) = \{\Y.fn\^)^ 



(2) 



i=l 



where i is the site index, M is the total number of sites, \n)i is the n-particle occupation- 
number state at site i, and /^ is the corresponding complex amplitude, with Ylm l/nP = 1- 
This ansatz does not capture the correlations found in analogs of fractional quantum Hall 
states, and hence we do not explore this model in that regime (which is characterized by 
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the interparticle spacing being of order the magnetic length). As a mean-field approach, 



36|, 

37|. 



it also does not capture quantum fiuctuation effects seen near the critical point [35 . 
nor is it appropriate for describing strongly-correlated macroscopic superposition states 
Furthermore it neglects all short-range correlations. Despite these limitations, the Gutzwiller 

1. It 



35 



36 



approach compares well with exact methods, and strong coupling expansions 

has been used extensively to understand experimental results [l|, [ij, [l5|, and is well suited 

for studying the vortex physics that we consider here. 

Using equation ([2]) as a variational ansatz, we minimize the energy with respect to the 
{/*}. We then extract the density pi = X^n^l/nP ^^^ ^^e condensate order parameter 
«i = (o-i) = J2n V^ {fn-i) fn ^^ ^ach sitc. The condensate density p^ = |(aj)p is equal to 
the superfluid density in this model, and is generally not equal to the density. 

We use an iterative algorithm to determine the {fn} which minimize the energy. We use 
a square region with L sites per side with hard boundary conditions. We find that we must 
take L much larger than the effective trap radius so that our solutions do not depend on 
those boundaries. Typically we use 40 < L < 90. We calculate {Hrbh) using equation 
(j2]). Minimizing {Hrbh) with respect to /" with the constraint ^„ |/^p — 1 = gives L^ 
nonlinear eigenvalue equations, one for each site. 



fc, nn of j 



{ak)Vmfi_iRjk + (aI)V"^ + l/m+i^fcj) + ( -^m 



U 



U 



^(^r) + -]m + Xj]fi = 0, 



(3) 



where the sum is over all nearest neighbors of site j, m is the particle-number index. A,- 



is a Lagrange multiplier, and Rjk = exp 



where i 



-1. We iteratively 



:/;Mr-A(r) 

solve these equations: first choosing a trial order-parameter field < o;^- k where aj = (aj); 
then updating it by aj^ = J2n V^fn-i ( \ '^j f ) /« ( i '^j f ) ' where p is the iteration 



index. Similar calculations were performed by Scarola and Das Sarma [30'] to analyze the 
case where the single-particle Mott state is surrounded by a rotating superfiuid ring. 

Since equation ([3]) is highly nonlinear, we find that the solution that this iterative al- 
gorithm converges to is sensitive to the initial state we use. This feature allows us to see 
the hysteretic effects described in the text. Experiments will see similar hysteresis, but the 
precise details will differ from our calculations (for example the critical frequencies for vortex 
entry and egress will be somewhat modified). 
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Detection 

To detect the vortices we propose time-of-flight imaging [l|, |33|], where at time t = one 
turns off the lattice and the harmonic trap, letting the cloud expand. After some fixed time 
t one then produces an absorption image of the cloud using a resonant laser beam. In a 
weakly interacting gas, the density profile is related to the momentum distribution of the 
gas. As we argue below, the vortices will be observable in the time-of-flight images. 

To illustrate, we present a simple model where we neglect interactions during the time 
of flight. As applied to atoms expanding from different sites, this approximation should be 
sound: by the time atoms from different sites overlap, the density is so low that they have 
negligable chance of scattering. On the other hand, atoms from the same site do scatter 
off one-another. These interactions have two effects: (1) atoms from sites with higher 
occupation will be moving faster (the interaction energy is converted into kinetic energy), 
and (2) the interactions introduce phase shifts which depend on atom number. Both of 
these effects will slightly degrade the contrast of the resulting interference pattern. When 
the interaction energies are small compared to the band spacing (a condition required for 
the single band Hubbard model to be applicable) these effects will be small. 

Within our approximation, calculating the expansion reduces to a series of single-particle 
problems. Taking the initial wavefunction to be given by equation ([2]), after time t the 
wavefunction will be 

M / [at(()l"\ 

i=l \ n I 

where aj(0) is the operator which annihilates the single-particle state in site i of the lattice. 



This operator evolves via the Heisenberg equation of motion, ihdtaj{t) = cijit), Hi 



Ij \oj, -Ofi-ce 



where 



if free is the Hamiltonian for non-interacting particles. This is equivalent to evolving the single 
particle state annihilated by aj(t) via the free Schrodinger equation. 

For this analysis we use the notation that r is a vector in the x — y plane, and z represents 
the coordinate in the perpendicular direction. We take the initial (Wannier) state at each 
site, (pi (r, z) , to be gaussian: 



[r,z) 
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2Ai 



(5) 



where A = -i/ , and A^ = 4/ with uJqsc and u± being the small oscillation frequencies 
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for each well. In the geometry we envision, oo\_ ^ oOosc- The wavefunctions at a time t after 
release of the trap are calculated by Fourier transforming 0j (r, z) to momentum space, then 
time evolving under Hi-ccc and finally Fourier transforming back to position space to arrive 
at 0i (r, 2;,t) = 0j (?, t) f{z,t), where the only thing we need to know about the transverse 
wavefunction f{z,t) is that it is normalized so J \f{z,t)\'^dz = 1. The in-plane wavefunction 
is 



(pi (r, t) 



X^ 



1/2 
2 I ^^P 



,7r (A^ + iht/m) 
The column density of the expanding cloud is then 

M 



r - r; 



2 (A2 + iht/m) 



(6) 



n r, 



t) = / (^(t) \^\r, z) ^(r, z) |^(t)) rf^; = ^ [n, - n,,J |0, (r, t)|2 + 



i=\ 



M 



^ ttj^i (r, t) 



i=l 



(7) 

where Ui {nc,i) is the number of atoms (condensed atoms) initially at site i, and ip{r, z) is 
the bosonic field operator annihilating an atom at position (?, z). 

In the long time limit where the expanded cloud is much larger than the initial cloud (ie. 
Dt = ht/mX ^ Rtf), this expression further simplifies, and one has 



n{r,t) = p{r,t)[{N-N,) + \A{r,t)\'] 
pir,t) = (^A^)-^-V^? 
A(r,t) = ^a,e--^/^'\ 



(8) 

(9) 

(10) 



where A^ and N^. are the total number of particles and condensed particles, respectively. 
The envelope, p{r,t), is a Gaussian, refiecting the Gaussian shape of the Wannier state. The 
incoherent contribution (A^ — Nc)p{r, t) 

has no additional structure. This is a consequence of the Gutzwiller approximation, which 
neglects short range correlations. Adding these correlations would modify the shape of the 
background, but it will remain smooth. 

The interference term has the structure of the envelope p{r^ t) multiplied by the modulus 
squared of the discrete Fourier transform of the superfiuid order parameter. The discrete 
Fourier transform can be constructed by taking the continuous Fourier transform of the 
product of a square array of delta-functions, and a smooth function which interpolates the 
superfiuid order parameter. The resulting convolution produces of a series of Bragg peaks, 
each of which have an identical internal structure which is the Fourier transform of the 
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interpolated superfluid order parameter. The vortices will be visible in the structure of 
these peaks. 

Vortices in real-space lead to vortices in reciprocal space. This result is clearest for "lowest 



Landau level" vortex lattices 39| which are expressible as an analytic function oi z = x + iy 
multiplied by a Gaussian of the form e^'^' ^^ , where w is a length scale which sets the cloud 
size. Aside from a scale factor and a rotation, the continuous Fourier transform of such a 
function is identical to the original. More generally, the topological charge associated with 
the total number of vortices is conserved in the expansion process. 

Figure [6] shows a sample expansion image in the long expansion limit, where up to a scaling 
of the spatial axes, the density pattern only depends on the {/^}'s from the Gutzwiller ansatz 
wavefunction, and the ratio between the size of the Wannier states and the lattice spacing 
(X/d). We use the initial state shown in fig. |U where t/U = 0.06. Using d = 410nm and 
hard axis lattice depth of SO-E^j, which are the experimental values in ref. [12], we find that 
Vo = 9.3Er, which gives A = 75 nm. 
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